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ABSTRACT 

We consider the long-standing problem of predicting the hierarchical clustering amplitudes 
Sp in the strongly non-linear regime of gravitational evolution. N-body results for the non- 
linear evolution of the bispectruni (the Fourier transform of the three-point density correlation 
function) suggest a physically motivated ansatz that yields the strongly non-linear behavior of 
the skewness, ^3, starting from leading-order perturbation theory. When generalized to higher- 
order [p > 3) polyspectra or correlation functions, this ansatz leads to a good description of 
non-linear amplitudes in the strongly non-linear regime for both scale-free and cold dark matter 
models. Furthermore, these results allow us to provide a general fitting formula for the non- 
linear evolution of the bispectrum that interpolates between the weakly and strongly non-linear 
regimes, analogous to previous expressions for the power spectrum. 

Subject headings: cosmology: theory, cosmology: large-scale structure of universe; methods: 
numerical; methods: analytical 
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Introduction 



Analytic understanding of gravitational clustering in the strongly non-linear regime remains an elusive 
goal in the study of large-scale structure. While the growth of structure in the weakly non-linear regime, 



where the density field contrast (5 < 1, is well described by cosmological perturbation theory (Juszkiewicz 



Bouchet fc Colombi 1993; Bernardeau 1994, 1994b; Baugh, Gaztahaga fc Efstathiou 1995; Gaztanaga fc 



Baugh 1995 ; Juszkiewicz et al. 1995 ; Lokas et al. 1995 ), corresponding progress has not been made in 



understanding clustering when S ^ 1. This lack of progress traces in part to the relative complexity and in- 
tractability of the coupled BBGKY hierarchy of equations for the phase space distribution functions (Peebles 
1980); moreover, with the development of high-resolution N-body simulations, the non- linear gravitational 
evolution can now be tracked directly. However, simulations provide only limited physical insight into the 
complex non-linear phenomena involved in gravitational clustering. To the extent possible, one would like 
a qualitative and ideally quantitative understanding of non-linear clustering from first principles. Analytic 
models can also provide a useful check on simulation results, which have their own intrinsic limitations and 
uncertainties. 

Recently, cosmological perturbation theory has been extended further into the non- linear regime through 
the inclusion of next-to-leading order corrections: one-loop perturbative calculations of the power spectrum 



(Makino et al. 1992; 


Lokas et al. 1996; [Scoccimarro & Frieman 1996 


), bispectrum (^coccimarro 1997|; 


Scoccimarro et al. 1998, hereafter SCFFHM) and skewness (3coccimarro & Frieman 1996a 




Scoccimarro 



1997; Fosalba & Gaztahaga 1998) have been found to be in excellent agreement with N-body results even in 
the regime where the variance = {S^{R)) ~ 1 for the bispectrum and as large as ^ 10 for the power 
spectrum (where S{R) is the density field contrast smoothed through a window of radius R). However, 
it is clear that perturbation theory cannot be pushed successfully beyond this point to calculate clustering 
amplitudes in the strongly non-linear regime. The perturbation series is at best asymptotic, and it is expected 
to break down on small scales, where ^ 1. More importantly, the single-stream fluid approximation, 
upon which the perturbative solutions are based, must fail on small scales where shell-crossing occurs, and 
be replaced with the BBGKY equations. 

Nevertheless, non-linear gravitational clustering does appear to display striking scaling properties, and 
some success has been achieved in combining perturbation theory with concepts such as stable clustering to 
make predictions about strongly non-linear behavior. The stable clustering hypothesis states that on very 
small scales, where clustering has reached virial equilibrium, the mean relative velocity between pairs should 
exactly cancel the Hubble expansion. When coupled with the equations of motion, this leads to general 
predictions for the growth of the hierarchy of p-point density correlation functions (Peebles 1980). For scale- 
free initial conditions, i.e., power spectrum Pi{k) ~ fc", in an Einstein-de Sitter universe (with f2m = 1), 
it follows from self-similarity that the slope of the non-linear two-point function, ^(r) ~ r~^ , is related to 
the spectral index n of the initial perturbations by 7 = 3(n + 3)/(n + 5) (Davis & Peebles 1977). For a 
range of initial spectra n, the results of N-body simulations are consistent with this self-similar solution in 
the strongly non-linear regime, ^ > 10 - 100 (Efstathiou et al. 1988; polombi, Bouchet, fc Hernquist 19*9^ , 
hereafter CBH; Jain 1997; Couchman fc Peebles 1998| ). In current models of structure formation, such as cold 
dark matter (CDM) and its variants, the initial power spectrum is not precisely scale-free. Nevertheless, the 
spectral index n ~ d\n P/ d In fc generally varies slowly enough with scale that one expects scaling arguments 
to provide useful insight into non-linear clustering. In this paper, we use such considerations to study the 
non-linear evolution of higher-order correlations. 



A classic problem in this context is the prediction of the strongly non-linear hierarchical amplitudes 
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Sp{R) = {SP{R)) / {S'^{R))P^^ . For scale-free initial conditions in an Einstein-de Sitter universe, the stable 
clustering hypothesis extended to higher-order correlations (Peebles 1980, Jain 1997) implies that the Sp 
should be constant, independent of R, in the strongly non-linear regime, ^ 1, but it does not say what 
the amplitudes are. 

This result motivated several papers which attempted to calculate the hierarchical amplitudes from the 
equations of motion. Observations of galaxy clustering at small scales and the self-similarity solution of the 
BBGKY equations motivate the so-called hierarchical model for the connected p-point correlation function 
(Fry 1984b), 

Kp{xi,...,Xp)= {6{xi),...,S{Xp))c^^Qp,a Y[^'^B=QpY^ Y1 n 

a=l labclings edges a—1 labelings edges 

The product is over p—1 edges that link p galaxies (vertices) A,B, with a two-point correlation func- 
tion S^xY assigned to each edge. These configurations can be associated with 'tree' graphs, called p-trees. 
Topologically distinct p-trees, denoted by a, in general have different amplitudes, denoted by Qp,a, but those 
configurations which differ only by permutations of the labels l,...,p (and therefore correspond to the same 
topology) have the same amplitude. There are tp distinct p-trees (^3 — 1, t4—2, etc., see Fry (1984b) and 
Boschan, Szapudi & Szalay 1994) and a total of pP~^ labeled trees. The hierarchical model represents the 
connected p-point functions as sums of products of {p — 1) two-point functions, introducing at each level 
only as many extra parameters Qp,a as there are distinct topologies. In what we shall call the degenerate 
hierarchical model, the amplitudes Qp,a are furthermore independent of scale and configuration. In this case, 
Qp,a = Qp, and the hierarchical amplitudes Sp ~ pP~'^ Qp. In the general case, the amplitudes Qp depend 
on overall scale and configuration. For example, for Gaussian initial conditions, in the weakly non-linear 
regime, ct^ ^ 1, perturbation theory predicts a clustering pattern that is hierarchical but not degenerate. 

Using the BBGKY hierarchy and assuming a hierarchical form similar to Eq. (^ for the phase-space 
p-point distribution function, in the stable clustering limit Fry (1982, 1984) obtained {p > 3) 



in this case, different tree diagrams all have the same amplitude, i.e., the clustering pattern is degenerate. 
On the other hand, Hamilton (1988), correcting an unjustified symmetry assumption in Fry (1982, 1984), 
instead found 



Qp, snake — , Qp, star — (3) 

where "star" graphs correspond to those tree graphs in which one vertex is connected to the other {p—1) 
vertices, the rest being "snake" graphs. Summed over the snake graphs, (0) yields 



p\ /0-,\p-^ 
p 

Unfortunately, as emphasized by Hamilton (1988), these results are not physically meaningful solutions to 
the BBGKY hierarchy, but rather a direct consequence of the assumed factorization in phase-space. As a 
result, this approach leads to unphysical predictions such as that cluster-cluster correlations are equal to 
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galaxy-galaxy correlations to all orders. It remains to be seen whether physically relevant solutions to the 
BBGKY hierarchy which satisfy Eq. (|l|) really do exist. Despite these shortcomings, the results in Eq. (^) 
and Eq. ^ are often quoted in the literature as physically relevant solutions to the BBGKY hierarchy! 

The complexity of the BBGKY equations has prompted a phenomenological approach to the description 
of correlation functions based on the hierarchical model. An example is the factorizable hierarchical model 
(Bernardeau & Schaeffer 1992), which is completely specified by the star amplitudes. In this case the snake 
amplitudes at any given order are determined by the lower order star amplitudes, Qp,a — ^ii'f^^'^K where 
are the vertex weights for i lines and di (a) is the number of such vertices in diagrams with topology denoted 
by a. This is analogous to the pattern that emerges from PT at large scales, although the parameters Qp^a 
are in this case taken to be constant, independent of scale and configuration, as in the spherical model. 

A more general framework than the hierarchical model is given by the scale-invariant model (Balian & 
Schaeffer 1989), in which the connected p-point function obeys the scaling law 



Kp{Xxi,...,XXp) = X ^'>'^ Kp{xi,...,Xp), (5) 

where 7 is the index of the two-point function, £,{r,t) ^ (r/t")~^, with a = 4/[3(n + 3)]. Eq. (||) is the 
self-similar solution to the BBGKY hierarchy, which holds in the case of stable clustering for scale-free 
initial conditions in an Einstein-de Sitter universe (Peebles 1980). The hierarchical model of Eq. (l|) satisfies 
Eq. (^, but the latter is more general and can be satisfied by other functional forms than Eq. (^). 

The best observational constraints on higher-order galaxy correlation functions in the non-linear regime 
currently come from angular surveys, although this situation will soon change dramatically with the advent of 
large redshift surveys such as the Two Degree Field (2dF) and Sloan Digital Sky Survey (SDSS). Measurement 
of the angular three-point function (Groth & Peebles 1977) and four-point function (Fry & Peebles 1978) in 
the Lick survey provided the first observational evidence for the scale- invariant model of Eq. (||). Moreover, 
these observations are consistent with the hierarchical model, with Qg, « 1.3 and Q4 « 3. More recently, 
Gaztahaga (1994) measured the Sp parameters (3 < p < 9) in the APM galaxy survey and found good 
agreement with the scale- invariant model at small scales. Szapudi & Szalay (1998) investigated the third- 
and fourth-order cumulant correlators in the APM and found Q3 « 1 and « 3, in good agreement with 
the Lick survey results. On the other hand, when decomposing the average four-point amplitude Q4 into the 
two different topologies, Szapudi & Szalay (1998) (see also Szapudi et al. 1995) found Ra — (94, snake = 3.7 
(for the snake topology) and Rb = (54, star — 0.8 (star topology), whereas Fry & Peebles (1978) obtained 
Ra — 2.5 ± 0.6 and Ri, — 4.3 ± 1.2 for the Lick survey. One must keep in mind, however, that these two 
measurements differ in the type of four-point configurations they considered (that is, they did not test the 
hierarchical ansatz in its full generality), so it remains possible that a finer measurement of the four-point 
function as a function of configuration could reconcile these discrepancies. Furthermore, in the EDSGC 
survey (which is based on a subset of the same photographic plates as used in the APM), Szapudi, Meiksin 
& Nichol (1996) found that Q3 — 2.0, Q4 = 7.3, substantially higher than in the APM (see Szapudi & 
Gaztanaga 1998 for a detailed comparison of higher order correlations in these two surveys). The upcoming 
2dF and SDSS redshift surveys will be able to probe higher-order correlations with unprecedented accuracy 
and should help resolve these issues ( |Colombi, Szapudi, fc Szalay 1998 ). We note, however, that testing 



the hierarchical model in these cases will require redshift distortions to be taken into account; in particular, 
the large internal velocity dispersion of galaxy clusters leads to a strong configuration-dependence of the Qp 
amplitudes at small scales ( Scoccimarro, Couchman &: Frieman 19981 ). In fact, this expected violation of the 



(degenerate) hierarchical ansatz is seen in the three-point function measured in the Las Campanas Redshift 
Survey (Jing & Borner 1998). 
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In this paper, we combine cosmological perturbation theory with the observed scahng behavior of the 
three-point function in N-body simulations to extract predictions for the non-hnear clustering amplitudes 
Sp for scale-free and CDM initial power spectra. The end result is a simple analytic expression for Sp which 
appears to be valid at ct^ > 10 when compared with N-body results. Along the way, we also provide fitting 
formulae for the non-linear evolution of the bispectrum, the Fourier transform of the three-point function, 
in the spirit of the two-point results of Hamilton et al. (1991), Jain, Mo & White (1995), and Peacock & 
Dodds (1994, 1996). 

In recent work, Colombi et al. (1997) take a very different approach to 'extending' perturbative ex- 
pressions for the Sp parameters into the non-linear regime. In leading order perturbation theory (hereafter 
PT), valid for cr^ <C 1, the Sp arc functions of the linear spectral index, n = —3 — dhiaf {R)/d\nR, where 
o-|(i?) denotes the linear variance at smoothing scale R (Bernardeau 1994). Colombi, et al. find that the 
non-linear evolution of the Sp can be fit with the same PT expressions provided they let the spectral index 
n appearing therein become a free parameter as a function of the non-linear variance, n neff(a'^). (Note 
that ricff is not the slope of the non-linear power spectrum.) The function r7.on(<T^), which depends on the 
initial spectrum, is extracted from, e.g., the measurement of S3{a^) in an N-body simulation, assuming that 
the perturbative expression for the skewness remains valid throughout the non-linear regime, i.e., by fitting 
the N-body results to S:i{a^) = (34/7) — (rieff + 3). Remarkably, they find that the extracted neff(CT^), when 
substituted into the PT expressions for the Sp for p > 3, provide a good fit to the non-linear (N-body) 
evolution of the higher-order moments as well. Similar behavior was noted for the EDSGC moments by 
Szapudi, Meiksin, and Nichol (1996). This extended perturbation theory (EPT) in principle provides a com- 
plete phenomenological description of one-point density field statistics, from linear to non-linear scales. (We 
describe it as a phenomenological model, because it is not a systematic development of perturbation theory 
for the non-linear regime; rather it is based upon the observation that the pattern of the Sp in the non-linear 
regime is related by a single parameter to the pattern in the weakly non-linear (perturbative) regime.) Note 
that EPT does not give a predictive prescription for the non-linear Sp for arbitrary initial conditions: if the 
initial conditions are changed, a new N-body simulation must be run in order to fit the function n{a^). In 
addition, it does not explictly describe the evolution of the multi-point functions. 

The model discussed below, which we have dubbed "hyperextended perturbation theory" (HEPT), 
also contains phenomenological elements: it does not provide a rigorous, first-principles calculation of the 
strongly non-linear clustering amplitudes. However, it is based upon a simple physical picture of the non- 
linear evolution of the A^— point functions that appears to hold quite generally. Thus, unlike extended 
perturbation theory, hyperextended perturbation theory yields an analytic prediction for the non-linear Sp 
for arbitrary initial conditions. 

This paper is organized as follows. Section 2 reviews the non- linear evolution of the bispectrum in 
PT and in N-body simulations and physically motivates the main ideas behind HEPT. In Section 3 we 
describe HEPT and present analytic results for the Sp parameters in the non-linear regime, comparing them 
to measurements in numerical simulations. A fitting formula for the non-linear evolution of the bispectrum 
from the weakly to the strongly non-linear regime is given in Section 4. Finally, Section 5 presents our 
conclusions. 



- 6 - 



2. Extracting the Essence of HEPT: Non-Linear Evolution of the Bispectrum 

We are interested in predicting the non-linear behavior of the higher-order correlations. For this purpose, 
we review here the non-linear evolution of the bispectrum, the Fourier transform of the three-point spatial 
correlation function, which is the lowest-order correlation function sensitive to phase information. 

It is useful to first recall the nomenclature of higher order correlations in Fourier space. Defining the 
Fourier transform of the density contrast field by 

the p'^ order polyspectrum. Bp, is defined by the expectation value 



{S{ki) . . . S{kp)) = [Sd]p Bp{ki, kp). (7) 

where = Soiki + . . . + kp). For p = 2, this defines the power spectrum, and for p = 3 the bispectrum. 
It is also convenient to define the Qp hierarchical amplitudes in Fourier space [see Eq. (|l|)] by {p > 2) 

Q ^ Bp{ki, . . . ,kp) 

J2a=l Slabclingsricdgcs^C^AB), 

where P{k) is the power spectrum, and the sum in the denominator is over all the tree diagrams, as in 
Eq. (|). For example, Q3 = B^I{PiP2 + P2P3 + Q4 = B^/lP^P^P^i + (3 permutations) + P1P12P4 + 

(11 permutations)], and so on. Here Pi = P{ki) and Pij = P{ki + kj). The hierarchical Sp{R) parameters 
at smoothing scale R are the one-point counterpart of the Qp amplitudes smoothed over a window of radius 
R: 



Sp{R) 



J Bp{ki, ...,kp)Wi...Wp [du]p d^h ■ ■ ■ <Pkp 



{5^[R))P-^ [J P{k)W{kR)^d^k]P-^ 

where Wi = W{kiR), with W{kR) the top-hat window function in Fourier space. 



(9) 



The non-linear evolution of Q3 displays the main features expected to hold for higher-order polyspectra 
as well. As its behavior has been rather thoroughly studied in both PT and N-body simulations, we shall use 
the three-point function as a model to extrapolate to higher-order polyspectra. The p > 3-order correlation 
functions become increasingly difficult to measure in numerical simulations for increasing p, so less is known 



about their non-linear evolution (however, see the studies of the four-point function by, e.g., Bromley 1994 , 
3uto & Matsubara 1994, Mimshi & Melott 1998). On the other hand, there is partial information from 



studies of the evolution of one-point statistics such as the Sp up to p = 10 (Baugh, Gaztaiiaga & Efstathiou 
1995| ; CBH). 



There are basically three qualitatively distinct regimes for the non- linear evolution of Q3: 

a) Tree-Level. At large scales, leading-order (tree-level) PT provides an excellent description of gravi- 
tational clustering. In this regime, the Qp amplitudes in Eq. (^) are independent of the overall amplitude 
of the power spectrum. For scale-free initial conditions, Pk{k) ~ fc", the PT Qp are scale-invariant but 
configuration-dependent, that is, they depend on the ratios ki/kj and the angles ki ■ kj/{kikj) (Fry 1984b). 
The resulting Sp parameters depend on scale only through the spectral index at the smoothing scale and its 



derivatives (Bernardeau 1994), therefore they are constants for scale-free initial conditions. 
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b) One-Loop. By carrying the perturbative expansion for the p— point functions to next-to-leading order 
in the density contrast field, one includes the one-loop corrections to the tree-level results. One-Loop PT 
describes the dependence of Qp and Sp on the amplitude of the power spectrum that appears at intermediate 
scales. Depending on the initial spectral index, one-loop corrections tend to enhance (n^ff > — 1-4) or 
reduce (ricff < — 1.4) the tree-level configuration dependence of Q3 (Scoccimarro 1997). For equilateral 
triangle configurations, one-loop PT describes the rise that connects the tree-level perturbative amplitude, 
Qs = 4/7, to the non- linear saturation value at small scales (SCFFHM). Similarly, one- loop corrections to 
the skewness 6*3 describe the transition from its tree- level value to the non-linear regime (Scoccimarro 1997). 
Extension of this result within the spherical collapse approximation shows that one-loop corrections describe 



this transition region for Sp with p > 3 as well ( Fosalba fc Gaztahaga 199 



c) Saturation. At small scales, in the strongly non-linear regime, numerical simulations show that the 
Sp{R) parameters reach a plateau, nearly independent of scale i?; we refer to this behavior as "saturation". 
The hierarchical amplitude in this regime becomes constant to a good approximation, nearly independent 
of configuration (SCFFHM). Note, however, that it is still not settled whether the Sp show small deviations 



from scale- invariant behavior in the strongly non- linear regime (CBH; Munshi et al. 1997). From a theoretical 
point of view, as we discussed above, the only prediction in the strongly non-linear regime is from stable 
clustering, which implies that the Sp should be constant for scale-free initial conditions in an Einstein-de 
Sitter model. As for the Qp parameters, however, stable clustering only constrains them to be scale-invariant 
(not necessarily hierarchical); in particular, for the bispectrum this implies 



S(fci, fc2, fcs) - P(fci)P(fc2) 5(ri2, 012) + P{k2)P{kz) 5(r23, 023) + P{kz)P[ki) 5(r3i, 03i). (10) 

Here, S{r,6) is some arbitrary function (symmetrized over ki and k2) of the ratios r^j = ki/kj and angles 
cosOij = {ki ■ kj)/{kikj); in the hierarchical model, S{r,9) = Q3 — constant. In order to preserve self- 
similarity, the time dependence of the bispectrum is driven by that of the power spectrum. The form in 



Eq. (10) generalizes to higher-order polyspectra, as in Eq. (|1|), where now each amplitude Qp^a becomes a 
different function of the ratios and angles 0y . This is exactly the behavior of higher-order correlations 
in tree-level PT (Fry 1984b). 

Since symmetries do not require a hierarchical three-point function in the non-linear regime, the N-body 
results suggest that the physics of gravitational clustering leads to S{r,9) being nearly constant, i.e., the 
saturation value for is only weakly dependent on r and 9. To see how this might arise, consider the 
interpretation of the Fourier hierarchical amplitude Q3 in tree-level PT. In this case, the function 5(r, 6) is 
given by 

10 . . / 1\ 4 2^ 



5(r,6i) = — -l-cos^ r+- +-cos^6i, (11) 
7 \ r / 7 

obtained from second-order PT (Fry 1984b). The configuration dependence through r and 9 in Eq. ( pi] ) 
comes from gradients of the density and velocity fields in the direction of the flow: the dependence on 
configuration arises from the anisotropy of structures and flows generated by the physics of gravitational 
instability (Scoccimarro 1997). Eq. ( pT|) implies that the hierarchical amplitude is maximum for coUinear 
configurations [cos 9 = ±1) and minimum for isosceles configurations (where two sides of the triangle are 
equal). This reflects the fact that gravitational instability generates large-scale flows mostly parallel to 
density gradients, which enhances coUinear configurations in Fourier space. On the other hand, on small 
scales, where virialization leads to substantial velocity dispersion, this picture suggests that the function 
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S{r, 9) should approach a constant. That is, non-colhnear configurations become more probable than at 
large scales: the loss of coherence between structures and flows implies that there is no reason to expect 
some configurations to be enhanced over others. 

Figure 1 illustrates these points for a CDM simulation, done by Couchman, Thomas & Pearce (1995) with 
an adaptive P'^M code that involves 128^ particles in a box of length 100 Mpc [h = i?o/100 kms^^ Mpc^^, 
where Hq is the Hubble constant). These simulation data are publicly available through the Hydra Consor- 
tium Web page ( http: / /coho.astro.uwo.ca/pub/consort.htm] ) and correspond to an il„i = 1 model, with linear 



CDM power spectrum characterized by a shape parameter F = 0.25, and normalization as = 0.64 (rCDM). 
The top right panel shows the dependence on scale of the hierarchical amplitude Qeq for equilateral triangles 
in Fourier space. At large scales (small k), Qeq approaches the tree- level PT value shown in the dashed 
line, Qeq = 4/7 = 0.57. At intermediate scales, Qeq rises {one-loop regime), eventually fiattening at small 
scales {saturation) . As shown in the figure, one-loop PT provides a good description of the transition regime 
but clearly breaks down in the saturation regime. The lower panels in Fig. 1 illustrate the dependence of 
the saturation value of Q3 on configuration, for fci/fc2 = 2 (bottom right) and fci/fc2 = 3,4 (bottom left). 
Although self-similarity considerations do not strictly apply to scale-dependent spectra such as the rCDM 
model, there is remarkably little dependence on configuration in the strongly non-linear regime (perhaps a 
slight decrease of Q3 with increasing fci/fc2 ratio), which confirms the expectations based on the physical 
picture discussed above. 

An important observation from Fig. 1 is that the saturation value of Q3 is in good agreement with the 
coUinear configuration value Q"'"'"(0, tt) given by tree- level PT (maxima of the dashed curves, at 6* = 0, tt). Al- 
though the dashed curves correspond to ^1/^2 = 2 configurations, the tree-level coUinear amplitude Q'^^(0, tt) 
is very insensitive to the ratio fci/fc2. In fact, for spectral indices n = —2,0, Q"^'"(0,7r) is independent of the 
ratio ki/k2- For these spectra, and in general to an excellent approximation, coUinear configurations are in- 
variant under arbitrary scaling of the different triangle sides. This property singles out these configurations, 
and we will take advantage of it to predict higher-order correlation amplitudes Sp in the non-linear regime. 



3. Hyperextended Perturbation Theory 

We shall now assume that clustering does reach approximate scale- invariance (saturation) at small scales 
(at least locally for CDM spectra). Based on the results above, we shall propose a physically motivated ansatz 
that allows one to calculate the Sp parameters in the non-linear regime purely from knowledge of tree-level 
PT. 

The discussion above shows that coUinear configurations play a special role in gravitational clustering. 
They correspond to matter flowing parallel to density gradients, thus enhancing clustering at small scales 
until eventually giving rise to bound objects that support themselves by velocity dispersion (virialization) . 
We thus conjecture that the "effective" Qp clustering amplitudes in the strongly non-linear regime are the 
same as the weakly non-linear (tree-level PT) coUinear amplitudes, as shown in Fig. 1 to hold for three-point 
correlations. We name this ansatz "hyperextended perturbation theory" (HEPT), since it borrows PT ideas 
valid at the largest scales to predict the behavior of clustering at small scales. 

Note that by effective amplitudes Qp^ we mean the overall magnitude of Qpi it is possible that Qp, for 
p > 3, although scale- invariant, is a function of configuration (as, e.g., in a non-degenerate hierarchical model, 
in which different topologies have different amplitudes Qp,a)- To calculate the resulting Sp parameters, we 
assume that Sp ~ pP-'^ Qf, that is, the Sp are given by the typical configuration amplitude Q^^ times the 
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total number of labeled trees, ^. In practice, there is a small correction to this formula due to smoothing, 
which we neglect ( [Boschan, Szapudi fc Szalay 1994 ). 



To obtain quantitative predictions, we must specify which tree-level collinear configurations we use to 
calculate. As p increases from 3, there is a growing number of collinear configurations, with different ratios 
Tij — ki/kj. However, as noted above, in tree-level PT different collinear amplitudes depend only very 
weakly on r^ , as in the p = 3 case. In what follows, we choose the fci = . . . = fcp_i = q, kp = —(j> — l)q 
configuration to calculate Qp. The resulting non- linear Sp amplitudes follow from tree-level PT 

4 — 2" 

sr(») = 3 ori") = 3 Y+¥^- 



-.sat/ N lor^sat. ^ .or 1536 - 1152 2" -f 128 3" -f 66 4" 64 6" - 9 8" - 2 12" - 24" 

? (Til 125 (D in] = 125 

^ ^ ' ^5 w 6 (1 + 12 2" + 12 3" + 16 4" + 24 6" -I- 24 8" -f 12 12" + 24 24") ' 



where n is the spectral index, obtained from (n + 3) = —d\nai{R)/d\nR, where cr^(i?) denotes the linear 
variance at smoothing scale R. One can check that these Qp amplitudes satisfy the constraint that cluster- 
cluster correlations are stronger than galaxy-galaxy correlations, 

^^^K^)'"^-^^--^2F^ (1^) 
( Hamilton fc Gott 1988| ), as long as n < 0.75, well within the physically interesting range. The constraint 
that the one-point probability distribution function is positive definite leads to Qp > pP"^ (Fry 1984), which 
is weaker than Eq. (|l^) and thus automatically satisfied. 

Figure 2 shows a comparison of these predictions with the numerical simulation measurements of CBH 
for scale-free initial conditions in an Einstein-de Sitter universe. These N-body simulations used a tree code 
(Hernquist, Bouchet fc Suto 1991) to evolve 64'^ particles in a cubic box with periodic boundary conditions. 
The plotted values correspond to the measured value of Sp when the non-linear variance cr^ — 100 (see CBH, 
table 3). The error bars denote the uncertainties due to the finite- volume correction applied to the raw 
Sp values. We see that the N-body results are generally in good agreement with the predictions of HEPT, 
Eqs. (|l2|), (13) and (^4|). The small discrepancy at n = —2 may be due to the excessive large-scale power in 



this model: in this case, the finite-volume corrections to the Sp measured in the simulations are quite large 
and thus uncertain. 



Figure 3 shows a similar comparison of HEPT with numerical simulations (Colombi, Bouchet, fc Scha- 



effer 1994 ) in the non-linear regime for the standard CDM model (with P = 0.5, as = 0.34). These N-body 
simulations used a P^M code to evolve 64^ particles in a box 64 Mpc on a side. The simulation error bars 
were estimated as in Fig. 2. The agreement between the N-body results and the HEPT predictions is very 
encouraging indeed. The change in the HEPT saturation value of the Sp with scale is due to the scale- 
dependence of the linear CDM spectral index, and follows the N-body results from ~ 10 to — 300, 
where stable clustering is approximately expected to hold. 

Is interesting to note that for n — 0, HEPT predicts Sp = {2p — 3)!!, which agrees exactly with the 
excursion set model developed by Sheth (1998) for white-noise Gaussian initial fluctuations. In this case, 
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the one-point probability distribution function for the density field yields an inverse Gaussian distribution, 
which has been shown to agree very well in the non-linear regime when compared to numerical simulations 
(Sheth 1998). This remarkable agreement between HEPT and the excursion set model deserves further work 
to understand the relation, if any, between these seemingly very different approaches to the description of 
statistics in the non-linear regime. 

From these comparisons, we conclude that HEPT provides a very good description of one-point statistics 
in the strongly non-linear regime. Note that this ansatz, although physically motivated, has not been 
rigorously proved from a theoretical point of view. Such a proof may be extremely difficult to achieve, 
due to the complexity of gravitational instability in the strongly nonlinear regime. However, the success of 
HEPT suggests that there is a deep connection between the physics of the non-linear regime and large-scale 
clustering. Such a connection was recently noted by Colombi et al. (1997) in formulating EPT, although 
the reason for it was not identified. 



4. A Fitting Formula for the Bispectrum 

In this section we provide a fitting formula that describes the non-linear evolution of the bispectrum as 
a function of scale, analogous to previous results in the literature for the power spectrum. The expression 
below interpolates between the perturbative results at tree-level (Fry 1984b) and one-loop (Scoccimarro 



1997; SCFFHM) and the saturation regime at small scales as studied by numerical simulations (Fry, Melott 



fc Shandarin 199^ ; SCFFHM) and described by HEPT. In order to describe both the weakly and strongly 



non-linear regimes, we take the following form for the bispectrum, inspired by the PT expression: 

B{ki,k2,k3) = 2Ff{ki,k2) P{ki) P(fc2) + permutations, (16) 

where P{k) is the non-linear power spectrum (obtained, e.g., from the fitting formulae of Peacock & Dodds 
1996), and the effective kernel 

Ff{ki,k2) = I a(7i, fci)a(n, fcz) + l ^i f ^ + ^) b{n,ki)b{n,k2) + ^ ^^^ c{n, ki)c{n, k2). (17) 
7 Z ki k2 ^k2 kiJ " ki K2 

When the functions a = 6 = c= l, we recover the tree-level PT expression for the bispectrum; on the other 
hand, for a? — (7/10) Ql*** and 5 = c = 0, we recover the results of HEPT for the strongly non-linear regime. 
The functions a(n, fc), b{n, fc), and c(n, k) are chosen to interpolate between these two regimes according to 
the one-loop PT and N-body resuhs (Scoccimarro 1997; SCFFHM). This yields (-2 < n < 0) 



1 + [0. 7 Qf^n)]^ /"^ {kRo)"+^ 
1 + [kRa 



= ^ + °f^: + ^i^^3f"" , (19) 
^ ' 1 + (fci?o)"+^-^ 



1 + 4.5/[L5jKn_f_3)4] (fci?o)"+^ 
1 + {kRoY 



= 'T/.rJ,:U. . (20) 



where the saturation value for the reduced bispectrum Q|^'(n) is given by Eq. (|T^), and Rq is the value of 
the correlation length in linear theory for Gaussian smoothing, i.e., cr|(i?o) = 1 with a Gaussian window 
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function. If desired, we can allow for residual scale-dependence of the saturation value by taking Ql 



Qsat (0-2^2OO)-'^'-"-' (Colombi et al. 1997; see also CBH). Here, we will focus on scale-free models and assume 



that scaling is achieved on the smallest scales. 

Figure 4 shows the fitting formula for the hierarchical three-point amplitude Qa (solid curves) as a 
function of angle 9 between fci and k2 for configurations with fci/fc2 = 2 configurations. These configurations 
are studied at different scales characterized by the value of kiR^. The power spectrum amplitude A(fc) = 
ATrk^P(k) provides a measure of the degree of non-linearity on these scales. The numerical simulation 
measurements shown here are taken from Figs. 1, 2, and 3 of SCFFHM. These N-body results correspond 
to 256^ PM simulations run by E. Hivon; the bispectrum measurements were done by S. Colombi, according 



to the scheme outlined in Appendix A of SCFFHM. In order to obtain from Eq. (16) we have used 
Eq. (|[) and used the linear power spectrum P{k) oc fc"; due to cancellations between the numerator and 
denominator, using the non- linear power spectrum (e.g., as given by Peacock fc Dodds 1996| ) does not change 



appreciably, though it does affect the bispectrum itself. Figure 5 shows a plot of the evolution of the 
equilateral hierarchical amplitude Qeq as a function of scale for n = —1.5 scale- free initial conditions. The 
results in Figs. 4 and 5 were found to be typical of the accuracy of the fitting formula, which is of order 15% 
for the scale-free models considered in SCFFHM. Generalization of Eqs. (p^)-(pO|) to CDM spectra is under 
way and will be reported elsewhere. 

Note that at intermediate scales, the functions a{n,k), b{n,k), and c{n,k) in Eq. ( pT|) break the scale 
invariance of Q, as required in the one-loop regime. We have chosen the overall scale dependence to be 
separable, in order to simplify applications of the fitting formula to large-scale structure calculations. For 
example, the skewness may be interpolated by 

53(i^)=f 4+H-4("+3), (21) 

7 cr* 7 cr* '^^ 



where 



aj{R) = / d'k P{k) j{n,k) W'{kR), (22) 



for j — a, b, c. In this way, only straightforward one-dimensional numerical integrations are required to 
calculate 6*3 at any scale R, using, e.g., the non-linear power spectrum from Peacock & Dodds (1996). 
Figure 6 shows an application of this result for n = — 1 scale- free initial conditions (solid curve), compared 
to numerical simulation measurements by CBH. Different symbols denote different expansion factors, a — 2.5 
(triangles), a = 6.4 (squares), and a = 16 (pentagons), where initial conditions were set at a = 1. Error bars, 
not shown for clarity, are of the order of 20% (CBH) ; thus the results of Eq. ( pl| ) are well within the N-body 
uncertainties. It seems, however, that the fit in Eq. ( ^ ) lies consistently above the N-body results. This is 
to be expected, at least up to cr^ « 1, since the numerical simulation measurements are affected by transients 
from the Zel'dovich approximation used to set up the initial conditions (Scoccimarro 1998). In fact, for the 
a = 2.5 output (triangles) in Fig. 6, the ct^ — > limit including transients corresponds to S3 = 2.39 instead 
of the tree-level PT value S3 = 2.86 shown by the dashed line. 



5. Conclusion 



We have proposed a simple ansatz, called Hyperextended Perturbation Theory (HEPT), to calculate 
clustering amplitudes, such as the hierarchical Sp parameters, in the strongly non-linear regime. Based on 
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N-body studies of scale-free and CDM models, the HEPT ansatz appears to be valid for all initial conditions 
of physical interest and has the advantage that it contains no free parameters. The proposal is based 
on extrapolating tree-level PT to the non-linear regime by using configurations that correspond to matter 
flowing parallel to density gradients. 

The similarity between the hierarchy of tree-level and strongly non-linear Sp parameters was pointed out 
by Colombi et al. (1997), who noticed that by using the tree- level expressions for the Sp with an arbitrary 
spectral index they could fit the whole hierarchy of Sp parameters in the non-linear regime. However, no 
explanation for such a remarkable coincidence was given. Within the HEPT framework, this arises naturally 
as a consequence of identifying the most important configurations that drive gravitational collapse. This 
provides a valuable tool for quantitative understanding of the behavior of higher-order correlations in the 
non-linear regime, based on simple physics borrowed from perturbation theory. At the same time, we 
caution that these results are only meant to provide insight into non-linear gravitational clustering: in the 
real universe, the clustering of luminous galaxies on these scales will be strongly affected by non-gravitational 
phenomena as well, such as gas dissipation, star formation, shocks, etc. 

We have also provided a fitting formula for the evolution of the bispectrum in scale-free models from 
the weakly to the strongly non-linear regime. This expression interpolates between the perturbative results 
on large scales and the saturation behavior observed in simulations at small scales. This result should be 
useful for applications in large-scale structure calculations. For example, using our results to obtain the 
skewness 6*3 as a function of scale and initial spectrum, one can implement the EPT framework of Colombi 
et al. (1997) to predict the rest of the one-point cumulants, Sp, for p > 3 at any scale. 

We thank Francis Bernardcau, Andrew Hamilton, and Istvan Szapudi for useful discussions and Stephane 
Colombi for providing the numerical simulation measurements in CBH and Colombi, Bouchet, & Schaeffer 
(1994). We also thank Ravi Sheth for pointing out the agreement of HEPT with the excursion set model 
for white-noise Gaussian fluctuations. This research was supported in part by the DOE at Chicago and 
Fermilab and by NASA grant NAG5-7092 at Fermilab. The tCDM simulations analyzed in this work 
were obtained from the data bank of cosmological A^-body simulations provided by the Hydra consortium 
://coho. astro. uwo.ca /pub/data. htm] ) and produced using the Hydra A^-body code (Couchman, Thomas, 
k Pearce 1995). 
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Fig. 1. — The upper left panel shows the power spectrum, A(fc) = 4Trk^P{k), for the rCDM model {as = 
0.64). Other panels show the hierarchical three-point Fourier amplitude Q3 for equilateral triangles (top 
right) and for configurations with different ki/k2 ratios (bottom) in the non-linear regime (in the lower left, 
triangles correspond to ki/k2 = 3, squares to ki/k2 = 4). Dashed curves show the predictions of tree-level 
PT. The solid curve in upper right panel shows the prediction of one-loop PT. 



- 16 - 




Fig. 2. — The hierarchical Sp parameters: skewness (triangles, p = 3), kurtosis (squares, p = 4), and pentosis 
(pentagons, p — 5) are shown for scale-free initial conditions in numerical simulations by Colombi, Bouchct 
& Hernquist (1996). The solid curves denote the predictions of HEPT, Eqs. (|l|), (H), and (|l|). 
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Fig. 3. — The skewness (triangles), kurtosis (squares), and pentosis (pentagons) parameters for the standard 
CDM spectrum (cts = 0.34) from numerical simulations by Colombi, Bouchet, & Schaeffer (1994). The solid 
curves denote the predictions of HEPT, Eqs. (|l|), (|l|) and (^. 
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Fig. 4. — The three-point hierarchical Fourier amphtude for configurations with ki/k2 = 2 as a function 
of the angle 9 between fci and k2, at different scales for n = —1.5 (top) and n = —2 (bottom) scale-free 



initial conditions. Symbols denote measurements in numerical simulations (taken from Scoccimarro et al 



1998, Figs. 1, 2, and 3). The solid curves show the bispectrum fitting formula, Eqs. ( pj[ ) to (|20|). The dashed 



curves correspond to the predictions of tree-level PT. 
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Fig. 5. — The three-point hierarchical Fourier amphtude Qeq for equilateral triangle configurations and 
n = —1.5 scale-free initial conditions as a function of scale. Symbols denote measurements in numerical 
simulations (taken from Scoccimarro et al. 1998 , Fig. 9). The solid curve shows the bispectrum fitting 
formula, Eqs. ( p^ ) to ([2C|), and the dashed line corresponds to the prediction of tree-level FT, Qeq = 4/7. 
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Fig. 6. — Predictions of our fitting formula for the skewness for n = — 1 scale- free initial conditions 
as a function of the variance for top-hat smoothing. Different symbols denote numerical simulation 
measurements at different expansion factors, a = 2.5 (triangles), a = 6.4 (squares), and a = 16 (pentagons), 
where initial conditions were set at a = 1. Error bars, not shown to avoid clutter, are of the order of 20% 
( polombi, Bouchet, fc Hernquist 1996 ). The solid curve shows the fitting formula of Eq. ([2]]), and the dashed 
line corresponds to the prediction of tree- level PT, S'3 = 20/7. The systematic underestimate of N-body 
results is likely due to transients from initial conditions (see text). 



